pacman::p_load(sf, tidyverse, tmap, sfdep, onemapsgapi, httr, jsonlite, olsrr, ggpubr, GWmodel, dotenv, matrixStats)Take-Home Exercise 3
Loading in all the packages and data
NOTE: TAKE THE PROF’s MPSZ 2019 DATA FROM INCLASS-EX 09!!!! Rather than 2014
Compile to rds once the data wrangling part is done.
Geographically weighted methods (REgression / Random FoREST) and Ordinary Least Square Method (Conventional non-geo-weighted linear regression)
Note that corrplot cannot be loaded (Therefore has to be loaded from package)
Geospatial Data:
Singapore National Boundary and Master Plan 2019 subzone!
Singapore National Boundary is a polygon feature data showing the national boundary!
Master Plan 2019 subzone are information on URA 2019!
mpsz <- st_read(dsn="data/geospatial", "MPSZ-2019")Reading layer `MPSZ-2019' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS: WGS 84national_boundary <- st_read(dsn="data/geospatial", "CostalOutline")Reading layer `CostalOutline' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21Aspatial Data
Resale Flat prices
Training dataset from 1st January 2021 to 31st December 2022
Test dataset from 1 January to the last day of February 2023 resale prices For our
For our case, we will be using four-room flats for our prediction
resale_flat <- read_csv("data/aspatial/resale-flat-prices-from-jan-2017-onwards.csv")Locational Factors:
For our assignment, we will need to look at location for
Proximity to CBD
Proximity to Eldercare
Proximity to Food court/Hawker
Proximity to MRT
Proximity to Park
Proximity to ‘Good’ primary schools
Proximity to Shopping Malls
Proximity to supermarket
Number of Kindergartens within 350m
Number of childcare centers within 350m
Number of primary schools within 1km
We start by sourcing some of this data
Geospatial Locational Source:
MRT/LRT and Bus Stop Locations
We can source from: LTA DataMall MRT/LRT Shapefile (It’s in Polygon!)
Bus Stop: LTA DataMall Bus Stop
mrt_lrt <- st_read(dsn="data/geospatial", layer="RapidTransitSystemStation")Reading layer `RapidTransitSystemStation' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 220 features and 4 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 6068.209 ymin: 27478.44 xmax: 45377.5 ymax: 47913.58
Projected CRS: SVY21bus_stop <- st_read(dsn="data/geospatial", layer="BusStop")Reading layer `BusStop' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 5159 features and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21There is something wrong, the MRT/LRT dataset is giving us polygon type rather than point. Despite the dataset saying it’s point.
OneMapSG Service API
We can extract the following according to their theme on OneMapSG API
Childcare (childcare)
Eldercare (eldercare)
Hawker Center (Queryname: hawkercentre)
Kindergarten (Queryname: Kindergartens)
Parks (Queryname: nationalparks)
Extra if we have time to think about:
Libraries (Queryname: libraries) [NUMBER WITHIN 350M]
Integrated Screening Programmes (queryname: moh_isp_clinics)
Tourist Attractions (Queryname: tourism) [NUMBER WITHIN 1KM]
Process of going through to create the shp file:
Courtesy of Megan’s work. (The following code chunk will be repeated to create the shp file for all the themes above).
#load_dot_env(file=".env")
#token <- Sys.getenv("TOKEN")
#themetibble <- get_theme(token, "themename")
#themetibble
#themesf <- st_as_sf(themetibble, coords=c("Lng", "Lat"), crs=4326)
#st_write(themesf, "themename.shp")Load in all the shp data
childcare <- st_read(dsn="data/geospatial", layer="childcare")Reading layer `childcare' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 1925 features and 5 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
Geodetic CRS: WGS 84eldercare <- st_read(dsn="data/geospatial", layer="eldercare")Reading layer `eldercare' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 133 features and 4 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.7119 ymin: 1.271472 xmax: 103.9561 ymax: 1.439561
Geodetic CRS: WGS 84hawker_centre <- st_read(dsn="data/geospatial", layer="hawkercentre")Reading layer `hawkercentre' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 125 features and 18 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.6974 ymin: 1.272716 xmax: 103.9882 ymax: 1.449017
Geodetic CRS: WGS 84kindergarten <- st_read(dsn="data/geospatial", layer="kindergartens")Reading layer `kindergartens' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 448 features and 5 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.6887 ymin: 1.247759 xmax: 103.9717 ymax: 1.455452
Geodetic CRS: WGS 84parks <- st_read(dsn="data/geospatial", layer="nationalparks")Reading layer `nationalparks' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 421 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.6929 ymin: 1.214491 xmax: 104.0538 ymax: 1.462094
Geodetic CRS: WGS 84libraries <- st_read(dsn="data/geospatial", layer="libraries")Reading layer `libraries' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 31 features and 13 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.7045 ymin: 1.263922 xmax: 103.9494 ymax: 1.448197
Geodetic CRS: WGS 84isp_clinics <- st_read(dsn="data/geospatial", layer="moh_isp_clinics")Reading layer `moh_isp_clinics' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 378 features and 15 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.6907 ymin: 1.26397 xmax: 103.9903 ymax: 1.456037
Geodetic CRS: WGS 84tourism <- st_read(dsn="data/geospatial", layer="tourism")Reading layer `tourism' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 109 features and 15 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.684 ymin: 1.223 xmax: 103.974 ymax: 1.447
Geodetic CRS: WGS 84External Sources
Supermarkets (Source: Dataportal)
Primary School (Wikipedia, Reverse Geocoded using OneMap API)
Top Primary Schools (We will pick the top 10 based off SchLah’s dataset)
Malls (Taken from this dataset) [NOTE: SVY21]
supermarkets <- st_read("data/geospatial/supermarkets.kml")Reading layer `SUPERMARKETS' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial\supermarkets.kml'
using driver `KML'
Simple feature collection with 526 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 103.6258 ymin: 1.24715 xmax: 104.0036 ymax: 1.461526
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84Load Primary School Data with Lat/Long
The lat/long are taken from OneMapSG API
Top 10 schools as well, dervied from schlah
primary_sch <- read_csv("data/aspatial/primary_school_geo.csv")
top10_pri_sch <- read_csv("data/aspatial/top_10_primary_school_geo.csv")Load Shopping Malls
shopping_mall <- read_csv("data/aspatial/shopping_mall.csv")Data Wrangling (Geospatial)
List of task to do:
Convert some remaining data to Geospatial
Convert multipoint to point data (Removing the ‘z’) for supermarket
Convert all if not already in SVY21 into SVY21 (3414)
Remove all unnecessary columns
Check for null values
Convert all of our datasets from CSV to Geospatial (shp)
Primary Schools
Shopping Malls
Within primary school, the OneMapAPI is not able to properly return the correct coordinate as the name is ‘CATHOLIC HIGH SCHOOL (PRIMARY)’. For it to return the correct coordinate, the name has to be truncated without (PRIMARY).
primary_sch_sf <- st_as_sf(primary_sch, coords=c("LONG", "LAT"), crs=4326)
top_primary_sch_sf <- st_as_sf(top10_pri_sch, coords=c("LONG", "LAT"), crs=4326)shopping_mall_sf <- st_as_sf(shopping_mall, coords=c("longitude", "latitude"), crs=4326)Remove ‘Z-Dimension’ for Supermarket Data
supermarkets <- st_zm(supermarkets)Check and change all EPSG for Geospatial data to ESPG 3414
childcare3414 <- st_transform(childcare, crs=3414)
eldercare3414 <- st_transform(eldercare, crs=3414)
hawker_centre3414 <- st_transform(hawker_centre, crs=3414)
kindergarten3414 <- st_transform(kindergarten, crs=3414)
parks3414 <- st_transform(parks, crs=3414)
libraries3414 <- st_transform(libraries, crs=3414)
isp_clinics3414 <- st_transform(isp_clinics, crs=3414)
tourism3414 <- st_transform(tourism, crs=3414)
primary_sch_sf_3414 <- st_transform(primary_sch_sf, crs=3414)
top_primary_sch_sf_3414 <- st_transform(top_primary_sch_sf, crs=3414)
shopping_mall_sf_3414 <- st_transform(shopping_mall_sf, crs=3414)Check others:
st_crs(primary_sch_sf_3414)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]st_crs(top_primary_sch_sf_3414)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]st_crs(shopping_mall_sf_3414)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]st_crs(supermarkets)Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]Looks like Supermarket is still in WGS84, let’s fix that
supermarkets3414 <- st_transform(supermarkets, crs=3414)Check to make sure…
st_crs(supermarkets3414)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]Remove all unnecessary columns
After glimpsing at the sf dataframe, we are going to remove all useless columns.
Information such as addresses/descriptions are not necessary, as all we need are the point data
We will keep the name for identification purposes
We will keep name and geometry
childcare3414 <- childcare3414 %>%
select(c('NAME', 'geometry'))Check result
childcare3414Simple feature collection with 1925 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 11810.03 ymin: 25596.33 xmax: 45404.24 ymax: 49300.88
Projected CRS: SVY21 / Singapore TM
First 10 features:
NAME geometry
1 APOLLO INTERNATIONAL PRESCHOOL PRIVATE LIMITED POINT (40985.94 33848.38)
2 APPLE TREE PLAYHOUSE PTE LTD POINT (28308.65 45530.47)
3 Appleland Montessori Child Care Centre Pte Ltd POINT (17828.84 36607.36)
4 APPLELAND PLAYHOUSE POINT (25579.73 29221.89)
5 APRICOT ACADEMY (LAGUNA) PTE. LTD. POINT (38981.02 32483.41)
6 Arise Preschool POINT (21588.47 36307)
7 Artemis Preskool @ Tampines Pte Ltd (CC) POINT (39239.78 37501.4)
8 Artemis Preskool @ Woodleigh POINT (32389.52 35403.72)
9 ARTS JUNIOR MONTESSORI LLP POINT (25554.36 30090.82)
10 Arts Kidz Pre-School Pte Ltd POINT (28004.17 28442)Same thing, name and geometry
eldercare3414 <- eldercare3414 %>%
select(c('NAME', 'geometry'))Check result
eldercare3414Simple feature collection with 133 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 14481.92 ymin: 28218.43 xmax: 41665.14 ymax: 46804.9
Projected CRS: SVY21 / Singapore TM
First 10 features:
NAME
1 Yuhua Senior Activity Centre
2 THK SAC @ Kaki Bukit
3 THK SAC @ Boon Lay
4 PEACE-Connect Senior Activity Centre@5
5 THK SAC @ Beo Crescent
6 Silver ACE @ Bukit Merah
7 Lions Befrienders Senior Activity Centre @ Tampines Blk 499C
8 Care Corner Senior Activity Centre (WL569)
9 Fei Yue Senior Activity Centre (Bukit Batok Branch)
10 COMNET Senior Activity Centre @ 182 Rivervale Crescent
geometry
1 POINT (16614.08 36639.12)
2 POINT (38803.81 35098.77)
3 POINT (14481.92 36357.61)
4 POINT (31505.35 31853.52)
5 POINT (27218.35 30135.49)
6 POINT (27278.94 29350.17)
7 POINT (41665.14 37956.92)
8 POINT (23147.94 45761.17)
9 POINT (18820.58 36396.32)
10 POINT (36446.37 41376.9)hawker_centre3414 <- hawker_centre3414 %>%
select(c('NAME', 'geometry'))hawker_centre3414Simple feature collection with 125 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 12874.19 ymin: 28355.97 xmax: 45241.4 ymax: 47850.43
Projected CRS: SVY21 / Singapore TM
First 10 features:
NAME geometry
1 Market Street Hawker Centre POINT (29874.82 29650.7)
2 Marsiling Mall Hawker Centre POINT (22042.51 46139.03)
3 Margaret Drive Hawker Centre POINT (24816.7 31094.91)
4 Fernvale Hawker Centre & Market POINT (32867.9 41500.77)
5 One Punggol Hawker Centre POINT (35955.52 43336.13)
6 Bukit Canberra Hawker Centre POINT (26794.39 47850.43)
7 Senja Hawker Centre POINT (19953.85 41008.06)
8 Buangkok Hawker Centre POINT (34575.37 40482.99)
9 Bukit Batok West Hawker Centre POINT (17817.9 37478.45)
10 Telok Blangah Hawker Centre & Market POINT (25191 28414.9)isp_clinics3414 <- isp_clinics3414 %>%
select(c('NAME', 'geometry'))isp_clinics3414Simple feature collection with 378 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 12130.5 ymin: 27388.9 xmax: 45475.65 ymax: 48626.7
Projected CRS: SVY21 / Singapore TM
First 10 features:
NAME geometry
1 HEALTHWAY MEDICAL CLINIC POINT (35625.67 40809.19)
2 HEALTHWAY MEDICAL CLINIC POINT (31090.4 33520.89)
3 HEALTHWAY MEDICAL CLINIC POINT (40903.43 39446.08)
4 HEALTHWAY MEDICAL CLINIC POINT (41384.47 37152.14)
5 HEALTHWAY MEDICAL CLINIC POINT (27371.14 45806.89)
6 HEALTHWAY MEDICAL CLINIC POINT (28014.68 45473.31)
7 HEALTHWAY MEDICAL CLINIC POINT (26462.45 47759.61)
8 HEALTHWISE MEDICAL CLINIC & SURGERY POINT (33725.38 32238.48)
9 HL CLINIC PTE LTD POINT (26451.91 28874.8)
10 HO MEDICAL CENTRE PTE LTD POINT (21103.55 32951.85)kindergarten3414 <- kindergarten3414 %>%
select(c('NAME', 'geometry'))kindergarten3414Simple feature collection with 448 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 11909.7 ymin: 25596.33 xmax: 43395.47 ymax: 48562.06
Projected CRS: SVY21 / Singapore TM
First 10 features:
NAME
1 PCF Sparkletots Preschool @ Cheng San-Seletar Blk 435 (KN)
2 PCF Sparkletots Preschool @ Cheng San-Seletar Blk 533 (KN)
3 PCF Sparkletots Preschool @ Cheng San-Seletar Blk 556 (DS)
4 PCF Sparkletots Preschool @ Chong Pang Blk 107 (KN)
5 PCF Sparkletots Preschool @ Chong Pang Blk 122 (KN)
6 PCF Sparkletots Preschool @ Chua Chu Kang Blk 10 (KN)
7 PCF Sparkletots Preschool @ Clementi Blk 330 (DS)
8 PCF Sparkletots Preschool @ Eunos Blk 616 (KN)
9 PCF Sparkletots Preschool @ Fengshan Blk 126 (DS)
10 PCF Sparkletots Preschool @ Fernvale Blk 416A (DS)
geometry
1 POINT (30325.45 38859.25)
2 POINT (30190.51 39574.18)
3 POINT (30705.05 39337.9)
4 POINT (27354.73 45992.92)
5 POINT (27755.87 46300.26)
6 POINT (19307.62 40271.08)
7 POINT (20706.39 32892.83)
8 POINT (37089.3 34892.34)
9 POINT (39752.13 34487.87)
10 POINT (33190.98 41392.37)libraries3414 <- libraries3414 %>%
select(c('NAME', 'geometry'))libraries3414Simple feature collection with 31 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 13665.24 ymin: 27383.57 xmax: 40922.89 ymax: 47759.75
Projected CRS: SVY21 / Singapore TM
First 10 features:
NAME geometry
1 Choa Chu Kang Public Library POINT (18175.64 40767.45)
2 Punggol Regional Library POINT (36014.52 43461.09)
3 Ang Mo Kio Public Library POINT (29364.62 39642.82)
4 Bedok Public Library POINT (38947.87 34357.48)
5 Bishan Public Library POINT (29727.41 36885.15)
6 Bukit Batok Public Library POINT (18646.83 36869.57)
7 Bukit Panjang Public Library POINT (20324.67 40211.98)
8 Central Public Library POINT (30327.79 31087.66)
9 Cheng San Public Library POINT (34728.1 39397.71)
10 Clementi Public Library POINT (20319.52 33038.5)We only keep STN_NAM_DE and geometry
Note that the file is in polygon format, rather than point. Might have to deal with it.
https://stackoverflow.com/questions/52522872/r-sf-package-centroid-within-polygon
mrt_lrt <- mrt_lrt %>%
select(c('STN_NAM_DE', 'geometry'))Check
mrt_lrtSimple feature collection with 220 features and 1 field
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 6068.209 ymin: 27478.44 xmax: 45377.5 ymax: 47913.58
Projected CRS: SVY21
First 10 features:
STN_NAM_DE geometry
1 ESPLANADE MRT STATION POLYGON ((30566.07 30621.21...
2 PAYA LEBAR MRT STATION POLYGON ((34495.6 33384.44,...
3 DHOBY GHAUT MRT STATION POLYGON ((29293.51 31312.53...
4 DAKOTA MRT STATION POLYGON ((34055.08 32290.62...
5 LAVENDER MRT STATION POLYGON ((31236.5 32085.76,...
6 RENJONG LRT STATION POLYGON ((34382.66 40949.64...
7 DOVER MRT STATION POLYGON ((21987.25 32576.91...
8 HOUGANG MRT STATION POLYGON ((34585.74 39132.47...
9 PHOENIX LRT STATION POLYGON ((19602.92 40048.64...
10 ALJUNIED MRT STATION POLYGON ((33422.35 33158.32...parks3414 <- parks3414 %>%
select(c('NAME', 'geometry'))parks3414Simple feature collection with 421 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 12374.75 ymin: 21917.81 xmax: 52533.09 ymax: 49296.46
Projected CRS: SVY21 / Singapore TM
First 10 features:
NAME geometry
1 ENG KONG PLACE FC POINT (20647.02 35214.35)
2 JANGGUS GARDEN POINT (28409.23 48791.68)
3 JLN LIMAU MANIS PG POINT (40880.22 34107.22)
4 GARDEN VIEW PG POINT (31600.01 38049.84)
5 THOMSON GREEN PG POINT (27919.96 40074.49)
6 JLN RIANG PG POINT (31789.37 36584.45)
7 MEI HWAN CRESCENT PG POINT (31162.69 37212.14)
8 FULTON AVE PG POINT (28026.6 38242.5)
9 MIMOSA TERRACE PG POINT (31016.83 40664.48)
10 JLN GENENG INTERIM PK POINT (33295.7 37270.28)primary_sch_sf_3414 <- primary_sch_sf_3414 %>%
select(c('Primary school', 'geometry'))primary_sch_sfSimple feature collection with 181 features and 5 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.6878 ymin: 1.274958 xmax: 103.9628 ymax: 1.456608
Geodetic CRS: WGS 84
# A tibble: 181 × 6
`Primary school` Type Gender Locat…¹ TAG geometry
* <chr> <chr> <chr> <chr> <chr> <POINT [°]>
1 Admiralty Primary School Gove… Mixed Woodla… ADMI… (103.8 1.442635)
2 Ahmad Ibrahim Primary S… Gove… Mixed Yishun AHMA… (103.8329 1.433153)
3 Ai Tong School Gove… Mixed Bishan AI T… (103.833 1.360583)
4 Alexandra Primary School Gove… Mixed Bukit … ALEX… (103.8244 1.291334)
5 Anchor Green Primary Sc… Gove… Mixed Sengka… ANCH… (103.8872 1.39037)
6 Anderson Primary School Gove… Mixed Ang Mo… ANDE… (103.8414 1.384264)
7 Ang Mo Kio Primary Scho… Gove… Mixed Ang Mo… ANG … (103.8396 1.369322)
8 Anglo-Chinese School (J… Gove… Boys Central ANGL… (103.841 1.30935)
9 Anglo-Chinese School (P… Gove… Boys Bukit … ANGL… (103.8356 1.318371)
10 Angsana Primary School Gove… Mixed Tampin… ANGS… (103.9518 1.348553)
# … with 171 more rows, and abbreviated variable name ¹Locationsupermarkets3414 <- supermarkets3414 %>%
select(c('Name', 'geometry'))supermarkets3414Simple feature collection with 526 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 4901.188 ymin: 25529.08 xmax: 46948.22 ymax: 49233.6
Projected CRS: SVY21 / Singapore TM
First 10 features:
Name geometry
1 kml_1 POINT (35561.22 42685.17)
2 kml_2 POINT (32184.01 32947.46)
3 kml_3 POINT (33903.48 39480.46)
4 kml_4 POINT (37083.82 35017.47)
5 kml_5 POINT (41320.3 37283.82)
6 kml_6 POINT (41384.47 37152.14)
7 kml_7 POINT (30186.63 38602.77)
8 kml_8 POINT (28380.83 38842.16)
9 kml_9 POINT (34383.76 37311.19)
10 kml_10 POINT (29010.23 45755.51)tourism3414 <- tourism3414 %>%
select(c('NAME', 'geometry'))tourism3414Simple feature collection with 109 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 11380.23 ymin: 22858.67 xmax: 43659.54 ymax: 47627.69
Projected CRS: SVY21 / Singapore TM
First 10 features:
NAME
1 National Gallery Singapore
2 Sultan Mosque (Masjid Sultan) Singapore
3 Sri Mariamman Temple: Hindu Temple in Singapore
4 Armenian Church in Singapore
5 CHIJMES Singapore
6 St Andrewâ??s Cathedral- Singapore Architecture Landmark
7 Kreta Ayer Square
8 Albert Mall Trishaw Park
9 Chinatown Food Street
10 Chinatown Heritage Centre, Singapore
geometry
1 POINT (30007.41 30267.17)
2 POINT (30877.14 31594.08)
3 POINT (29342.34 29382.57)
4 POINT (29818.66 30598.9)
5 POINT (30043.47 30820.05)
6 POINT (30113.58 30488.32)
7 POINT (29142.6 29293.89)
8 POINT (30244.32 31323.04)
9 POINT (29165.65 29443.75)
10 POINT (29227.71 29603.72)Dealing with the anomaly - MRT/LRT
Despite it saying that’s point data, it gives polygon data, therefore we need to convert the polygon into point data. Unfortunately, there appears to be an issue that causes st_make_valid to not work, due to missing data within the polygon (There is an ‘na’).
Thus, we’ll take this chance to do some fine tuning and make use of OneMapSG to return us the correct coordinates. We will also remove all ‘depots’ as they are not considered MRT/LRTs for our dataset.
mrt_lrt_new_set <- read_csv("data/aspatial/mrt_lrt_lat_long.csv")Convert into an SF Object
mrt_lrt_ns_sf <- st_as_sf(mrt_lrt_new_set, coords=c("LONG", "LAT"), crs=4326)mrt_lrt_ns_sf <- mrt_lrt_ns_sf %>%
select(c('STN_NAM_DE', 'geometry'))mrt_lrt_ns_sf_3414 <- st_transform(mrt_lrt_ns_sf, crs=3414)View the map to make sure the mapping is correct
tmap_mode("view")+
tm_shape(bus_stop)+
tm_dots(col="purple", size=0.05)Looks like we have some bus stations in Malaysia, which makes sense, as there are some Singapore buses that travel into Malaysia.
These stops are:
LARKIN TERMINAL (Malaysia) [46239]
KOTARAYA II TER [46609]
JB SENTRAL [47701]
JOHOR BAHRU CHECKPOINT (46211, 46219)
We don’t want them in our calculation, so we’ll remove them.
bus_stop_to_remove <- c(46239, 46609, 47701, 46211, 46219)
bus_stop_cleaned <- filter(bus_stop, !(BUS_STOP_N %in% bus_stop_to_remove))Check the map again:s
tm_shape(bus_stop_cleaned)+
tm_dots(col="purple", size=0.05)Convert Bus Stop from WGS84 to SVY21
bus_stop_cleaned_3414 <- st_transform(bus_stop_cleaned, crs = 3414)st_crs(bus_stop_cleaned_3414)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]Data Wrangling (Aspatial)
We will strictly be looking at Four Room Flats
Training Set: 1st January 2021 - 31st December 2022
Test Set: January - February 2023
resale_flat_sub <- filter(resale_flat,flat_type == "4 ROOM") %>%
filter(month >= "2021-01" & month <= "2023-02")Transforming the Data into usable data
Looking at the data, we will now need to combine both the block and street name to form the address so we can lookup for the postal code and retrieve the LAT/LONG
We will also need to calculate the remaining number of lease. Given that the remaining lease are in ‘YEAR’ + ‘MONTH’ format, we need to change it into a continuous number. We will change it to ‘MONTHS’.
rs_modified <- resale_flat_sub %>%
mutate(resale_flat_sub, address = paste(block,street_name)) %>%
mutate(resale_flat_sub, remaining_lease_yr = as.integer(str_sub(remaining_lease, 0, 2))) %>%
mutate(resale_flat_sub, remaining_lease_mth = as.integer(str_sub(remaining_lease, 9, 11)))Calculate the LEASE in Months and remove all intermediary columns
Convert all NA into 0
Sum up YEAR and MONTH and place them into a new column
Select only the columns we need.
rs_modified$remaining_lease_mth[is.na(rs_modified$remaining_lease_mth)] <- 0
rs_modified$remaining_lease_yr <- rs_modified$remaining_lease_yr * 12
rs_modified <- rs_modified %>%
mutate(remaining_lease_mths = rowSums(rs_modified[, c("remaining_lease_yr", "remaining_lease_mth")])) %>%
select(month, town, address, block, street_name, flat_type, storey_range, floor_area_sqm, flat_model,
lease_commence_date, remaining_lease_mths, resale_price)Glimpse at our data:
glimpse(rs_modified)Retrieving the LAT/LONG Coordinates from OneMap.sg API
The following solution is courtesy of: Nor Aisyah
We first create a unique list of addresses
add_list <- sort(unique(rs_modified$address))Now we will focus on getting the data needed.
In the following code chunk, the goal is to retrieve coordinates for a given search value using OneMap SG’s REST APIs. The code creates a dataframe to store the final retrieved coordinates and uses the httr package’s GET() function to make a GET request to https://developers.onemap.sg/commonapi/search.
The required variables to be included in the GET request are searchVal, returnGeom {Y/N}, and getAddrDetails {Y/N}.
The JSON response returned will contain multiple fields, but the code is only interested in the postal code and coordinates like Latitude & Longitude. A new dataframe is created to store each final set of coordinates retrieved during the loop, and the number of responses returned is checked to append to the main dataframe accordingly. Invalid addresses are also checked, and the response is appended to the main dataframe using rbind() function of base R package.
retrieve_coords <- function(add_list){
# Create a data frame to store all retrieved coordinates
postal_coords <- data.frame()
for (i in add_list){
#print(i)
r <- GET('https://developers.onemap.sg/commonapi/search?',
query=list(searchVal=i,
returnGeom='Y',
getAddrDetails='Y'))
data <- fromJSON(rawToChar(r$content))
found <- data$found
res <- data$results
# Create a new data frame for each address
new_row <- data.frame()
# If single result, append
if (found == 1){
postal <- res$POSTAL
lat <- res$LATITUDE
lng <- res$LONGITUDE
new_row <- data.frame(address= i, postal = postal, latitude = lat, longitude = lng)
}
# If multiple results, drop NIL and append top 1
else if (found > 1){
# Remove those with NIL as postal
res_sub <- res[res$POSTAL != "NIL", ]
# Set as NA first if no Postal
if (nrow(res_sub) == 0) {
new_row <- data.frame(address= i, postal = NA, latitude = NA, longitude = NA)
}
else{
top1 <- head(res_sub, n = 1)
postal <- top1$POSTAL
lat <- top1$LATITUDE
lng <- top1$LONGITUDE
new_row <- data.frame(address= i, postal = postal, latitude = lat, longitude = lng)
}
}
else {
new_row <- data.frame(address= i, postal = NA, latitude = NA, longitude = NA)
}
# Add the row
postal_coords <- rbind(postal_coords, new_row)
}
return(postal_coords)
}We call the function to get the coordinates
coords <- retrieve_coords(add_list)Inspect the results
coords[(is.na(coords$postal) | is.na(coords$latitude) | is.na(coords$longitude) | coords$postal=="NIL"), ]From the results, we can see that it has missing postal code. However, the postal code is not as relevant to our analysis as latitude and longitude.
The missing respective postal code are:
680215
680216
rs_coords <- left_join(rs_modified, coords, by = c('address' = 'address'))Write file to RDS
write_rds(rs_coords, "data/aspatial/rs_coords.rds")Read the RDS file
resale_sub_flat <- read_rds("data/aspatial/rs_coords.rds")Assign and transform CRS
Data from OneMapSG are in WGS84, as evident from the decimal. So we will create a SF object in that code before transforming into SVY21
resale_coords_sf <- st_as_sf(resale_sub_flat,
coords = c("longitude",
"latitude"),
crs=4326) %>%
st_transform(crs = 3414)st_crs(resale_coords_sf)Check if there are any NA outside of what we confirmed earlier
rows_with_na <- which(is.na(resale_coords_sf), arr.ind=TRUE)if (nrow(rows_with_na) > 0) {
message("The following rows have NA values:")
print(rows_with_na)
} else {
message("The dataframe does not contain any rows with NA values.")
}Remove all unnecessary rows
We only need to keep the identifier, time, resale price, spatial point data, and the structural factors as listed below
Area of Unit
Floor Level
Remaining Lease
Age of Unit
trim_resale_flat <- resale_coords_sf %>%
select(1:3, 6:8, 11:12)Determining CBD
We need to consider the distance to CBD as well. We will take ‘Downtown Core’ as our reference point, which is located in the southwest of Singapore.
cbd_lat <- 1.287953
cbd_long <- 103.851784
cbd_sf <- data.frame(cbd_lat, cbd_long) %>%
st_as_sf(coords = c("cbd_long", "cbd_lat"), crs=4326) %>%
st_transform(crs=3414)Proximity Distance Calculation
For our next step, we will need to integrate all of our geospatial data by calculating the proximity to the amenities.
The following function calculates proximity
- rowMins from matrixStats package finds the shortest possible distance within the distance matrix.
proximity <- function(df1, df2, varname) {
dist_matrix <- st_distance(df1, df2) %>%
units::drop_units()
df1[,varname] <- rowMins(dist_matrix)
return(df1)
}Implement the Proximity Calculations:
trim_resale_flat <- proximity(trim_resale_flat, cbd_sf, "PROX_CBD") %>%
proximity(., eldercare3414, "PROX_ELDERCARE") %>%
proximity(., hawker_centre3414, "PROX_HAWKER") %>%
proximity(., isp_clinics3414, "PROX_ISP") %>%
proximity(., mrt_lrt_ns_sf_3414, "PROX_MRT") %>%
proximity(., parks3414, "PROX_PARKS") %>%
proximity(., top_primary_sch_sf_3414, "PROX_TOP_PRI") %>%
proximity(., shopping_mall_sf_3414, "PROX_SHOPPING") %>%
proximity(., supermarkets3414, "PROX_SUPERMARKETS")Facility Count within Radius
In addition to calculating the shortest distance between points, we’re also interested in finding out how many facilities are located within a specific radius. To accomplish this, we’ll use the st_distance() function to calculate the distance between the flats and the desired facilities. We’ll then use rowSums() to add up the observations and obtain the count of facilities within the desired radius. The resulting values will be added to the data frame as a new column.
num_radius <- function(df1, df2, varname, radius) {
dist_matrix <- st_distance(df1, df2) %>%
units::drop_units() %>%
as.data.frame()
df1[,varname] <- rowSums(dist_matrix <= radius)
return(df1)
}Implement the function:
trim_resale_flat <-
num_radius(trim_resale_flat, kindergarten3414, "NUM_KINDERGARTEN", 350) %>%
num_radius(., childcare3414, "NUM_CHILDCARE", 350) %>%
num_radius(., bus_stop_cleaned_3414, "NUM_BUS_STOPS", 350) %>%
num_radius(., primary_sch_sf_3414, "NUM_PRI_SCHS", 1000) %>%
num_radius(., tourism3414, "NUM_TOURIST_SITES", 1000)Save the Data as RDS
Before we do that, we should trim the name, as when converting to RDS the column names are shortened due to limitations
trim_resale_flat <- trim_resale_flat %>%
mutate() %>%
rename("AREA_SQM" = "floor_area_sqm",
"LEASE_MTHS" = "remaining_lease_mths",
"PRICE"= "resale_price",
"FLAT" = "flat_type",
"STOREY" = "storey_range")st_write(trim_resale_flat, "data/geospatial/resale_flat_final.shp")Reimport by reloading our dataset
resale_flat_sf <- st_read(dsn="data/geospatial", layer="resale_flat_final") Reading layer `resale_flat_final' from data source
`D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 25503 features and 22 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 11519.79 ymin: 28217.39 xmax: 42645.18 ymax: 48741.06
Projected CRS: SVY21 / Singapore TMWe’re done with gathering and creating our data! Now we can move onto our exploratory data analysis and building our prediction models. Note that there are still some data that need to be modified. Such as ordinal data like Storeys.
Speaking of which, we’ll handle them now.
Storey Data:
It’s a categorical data that has meaning when ordered. Higher/lower levels could potentially have an impact on the price of the HDB Flat. Therefore, we should convert this column of data into ordinal scale, one that the regression can learn from.